The following R packages are used in the analysis.
library(tidyverse)
library(DataExplorer)
library(brms)
library(rstanarm)
library(ggdist)
library(skimr)
library(Amelia)
library(caret)
library(rpart.plot)
library(gridExtra)
library(lme4)
library(BAS)
library(latex2exp)
# Set theme for ggplot
mytheme = theme_bw()
theme_set(mytheme) Read the data
d <- read_csv("PerturbedCyclistExperimentalDataTUDelft2022.csv",
col_types = cols(participant_ID = col_character(),
outcome = col_factor()))We reorder participant_ID according to age,
filter out initial search data and remove some columns that we don’t
further need in the analysis.
d <- d %>%
mutate(participant_ID = fct_reorder(participant_ID, age)) %>%
filter(initial_search==FALSE, exclude=="NO") %>%
select(-c(pull_force, comments, initial_search, exclude))
d %>% glimpse()## Rows: 928
## Columns: 19
## $ participant_ID <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ gender <chr> "female", "female", "female", "female", "female", "f…
## $ age <dbl> 61, 61, 61, 61, 61, 61, 61, 61, 61, 61, 61, 61, 61, …
## $ weight <dbl> 83.6, 83.6, 83.6, 83.6, 83.6, 83.6, 83.6, 83.6, 83.6…
## $ length <dbl> 1.61, 1.61, 1.61, 1.61, 1.61, 1.61, 1.61, 1.61, 1.61…
## $ skill_performance <dbl> 81.30124, 81.30124, 81.30124, 81.30124, 81.30124, 81…
## $ skill_effort <dbl> 3.798331, 3.798331, 3.798331, 3.798331, 3.798331, 3.…
## $ velocity <dbl> 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, …
## $ pull_ID <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1…
## $ pull_direction <chr> "CCW", "CCW", "CW", "CW", "CCW", "CCW", "CW", "CCW",…
## $ outcome <fct> 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 0…
## $ angular_momentum <dbl> 6.7600990, 9.2634520, 11.9062300, 2.9413360, 8.27263…
## $ reaction_time <dbl> 0.0045, 0.0190, 0.0760, 0.0990, 0.0490, 0.0655, 0.03…
## $ phi_zero <dbl> -1.176674000, -0.978592300, -1.416520000, -0.9767104…
## $ phi_dot_zero <dbl> 0.194354000, 0.163115900, -0.013537920, 0.080976880,…
## $ delta_zero <dbl> -0.030918070, -0.017670920, -0.018878270, -0.0093245…
## $ delta_dot_zero <dbl> 0.202679100, -0.127448200, -0.195207200, 0.156655100…
## $ psi_zero <dbl> -0.027438300, -0.003981016, 0.023122360, -0.01823309…
## $ Q_zero <dbl> -46.488590, -39.377390, 34.093770, 31.780850, -1.862…
d %>% skim()| Name | Piped data |
| Number of rows | 928 |
| Number of columns | 19 |
| _______________________ | |
| Column type frequency: | |
| character | 2 |
| factor | 2 |
| numeric | 15 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| gender | 0 | 1 | 4 | 6 | 0 | 2 | 0 |
| pull_direction | 4 | 1 | 2 | 3 | 0 | 2 | 0 |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| participant_ID | 0 | 1 | FALSE | 24 | 22: 60, 2: 60, 15: 60, 11: 60 |
| outcome | 0 | 1 | FALSE | 2 | 0: 484, 1: 444 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| age | 0 | 1.00 | 43.74 | 17.67 | 20.00 | 29.00 | 34.00 | 62.00 | 75.00 | ▇▃▁▅▃ |
| weight | 0 | 1.00 | 77.01 | 11.70 | 58.30 | 67.90 | 76.40 | 85.10 | 103.30 | ▆▇▇▅▂ |
| length | 0 | 1.00 | 1.81 | 0.09 | 1.61 | 1.73 | 1.83 | 1.87 | 1.98 | ▃▃▅▇▃ |
| skill_performance | 20 | 0.98 | 65.45 | 26.69 | 22.32 | 42.91 | 63.77 | 85.70 | 135.02 | ▇▇▇▅▁ |
| skill_effort | 0 | 1.00 | 4.78 | 3.56 | 0.00 | 2.57 | 4.18 | 6.42 | 18.82 | ▇▇▂▁▁ |
| velocity | 0 | 1.00 | 12.12 | 4.22 | 6.00 | 12.00 | 12.00 | 18.00 | 18.00 | ▃▁▇▁▅ |
| pull_ID | 0 | 1.00 | 22.49 | 14.72 | 1.00 | 10.00 | 20.00 | 33.00 | 60.00 | ▇▇▆▃▂ |
| angular_momentum | 0 | 1.00 | 13.95 | 5.22 | 0.00 | 10.37 | 13.80 | 17.54 | 29.23 | ▁▆▇▅▁ |
| reaction_time | 17 | 0.98 | 0.03 | 0.02 | 0.00 | 0.01 | 0.03 | 0.05 | 0.15 | ▇▇▂▁▁ |
| phi_zero | 24 | 0.97 | -0.98 | 1.14 | -4.41 | -1.85 | -0.92 | -0.13 | 4.10 | ▁▇▇▁▁ |
| phi_dot_zero | 15 | 0.98 | 0.03 | 0.14 | -0.92 | -0.01 | 0.03 | 0.07 | 1.04 | ▁▁▇▁▁ |
| delta_zero | 24 | 0.97 | 0.00 | 0.03 | -0.15 | -0.02 | 0.00 | 0.02 | 0.17 | ▁▂▇▁▁ |
| delta_dot_zero | 15 | 0.98 | 0.07 | 0.41 | -2.37 | -0.11 | 0.07 | 0.25 | 2.37 | ▁▁▇▁▁ |
| psi_zero | 24 | 0.97 | 0.00 | 0.02 | -0.09 | -0.01 | 0.00 | 0.01 | 0.12 | ▁▅▇▁▁ |
| Q_zero | 15 | 0.98 | 13.38 | 64.76 | -320.51 | -22.80 | 11.12 | 47.29 | 240.67 | ▁▁▇▆▁ |
plot_missing(d, missing_only=TRUE, ggtheme =mytheme)d %>% missmap(y.labels=d$participant_ID,
main="Missingness map, (participant_ID on vertical axis)")dcontinuous <- d %>% select(age, angular_momentum, delta_dot_zero, delta_zero, length,
phi_dot_zero, phi_zero, psi_zero, Q_zero, reaction_time,
skill_effort, skill_performance, weight)
dcontinuous %>% plot_histogram(ggtheme=mytheme)d %>% drop_na() %>%
mutate(`log10(reaction_time)`=log10(reaction_time)) %>%
ggplot(aes(x=`log10(reaction_time)`)) +
geom_histogram(bins=50,colour="white")plot_correlation(dcontinuous, ggtheme = mytheme)
Hence, only
length and weight have moderate
correlation.
fracfalling <- d %>% mutate(velocity_=as.factor(velocity)) %>%
group_by(participant_ID, velocity_) %>%
summarise(n=n(),
fraction_falling=mean(outcome=="1"),
skill_performance=mean(skill_performance),
participant_ID = first(participant_ID),
age=mean(age),
gender=first(gender),
skill_effort=first(skill_effort),
average_reaction_time=mean(reaction_time))
fracfalling## # A tibble: 47 × 9
## # Groups: participant_ID [24]
## participant_ID velocity_ n fraction_falling skill_performan… age gender
## <fct> <fct> <int> <dbl> <dbl> <dbl> <chr>
## 1 9 6 20 0.4 34.3 20 female
## 2 9 12 20 0.5 44.7 20 female
## 3 16 12 20 0.45 60.4 24 female
## 4 16 18 19 0.579 36.0 24 female
## 5 19 12 20 0.45 115. 27 female
## 6 19 18 20 0.45 86.3 27 female
## 7 22 6 20 0.35 46.3 27 male
## 8 22 12 20 0.5 75.0 27 male
## 9 22 18 20 0.5 58.0 27 male
## 10 6 12 20 0.55 85.7 27 male
## # … with 37 more rows, and 2 more variables: skill_effort <dbl>,
## # average_reaction_time <dbl>
Visualise
fracfalling %>%
ggplot(aes(x=participant_ID, y = fraction_falling, colour=velocity_)) +
geom_point() d %>%
ggplot(aes(x=angular_momentum, y=outcome)) +
geom_jitter(height=0.2)We can make separate panels for each participant and colour according to velocity
d %>% mutate(velocity=as.factor(velocity)) %>%
ggplot(aes(x=angular_momentum, y = outcome, colour=velocity)) +
geom_jitter(height=0.1,size=0.5) +
facet_wrap(~participant_ID,nrow=4) +
labs(x="angular momentum", y="y")+
theme(legend.position = "top")We investigate in which sense the variables phi_zero,
phi_dot_zero, delta_zero,
delta_dot_zero, psi_zero and
Q_zero differ over the binary outcome.
dmeas <- d %>% select(outcome, phi_zero, phi_dot_zero,
delta_zero, delta_dot_zero, psi_zero, Q_zero) %>%
pivot_longer(cols=contains("zero"), names_to="measurement", values_to="value")
dmeas %>%
ggplot(aes(sample=value,colour=outcome)) +
stat_qq() +
stat_qq_line() +
facet_wrap(~measurement, scales="free") +
labs(x="theoretical", y="sample")We observe similar patterns for all of these 6 variables.
d %>% ggplot(aes(sample=phi_dot_zero,colour=outcome, shape=pull_direction)) +
stat_qq(size=2)We fit a classification tree to get a first impression of variables that matter for predicting the outcome variable. The advantage of such a method is that it is insensitive to noise variables that don’t have predictive power.
tr <- train(outcome ~ ., data = na.omit(d), method='rpart', tuneLength=10)
rpart.plot(tr$finalModel,extra=101)ggplot(varImp(tr))Classification trees are highly variable. Therefore, we fit a random forest. As the fraction of rows containing missing values is relatively small, we fit the RF without rows containing missing values.
d_dropna <- d %>% drop_na()
tr_rf <- train(outcome ~ ., data = d_dropna, method='rf',
trControl=trainControl("cv"), importance = TRUE, tuneLength=10)
ggplot(varImp(tr_rf))We train the random forest again, but now such that participant_ID is not automatically converted to dummy variables. This implies that variable importance is computed for participant_ID by itself, and not for individual participants.
X <- d_dropna %>% dplyr::select(-outcome)
Y <- d_dropna %>% dplyr::select(outcome) %>% deframe()
tr_rf2 <- train(X,Y, method='rf',
trControl=trainControl("cv"), importance = TRUE, tuneLength=10)
p2 <- ggplot(varImp(tr_rf2))Refit, but now further excluding participant_ID.
Xsub <- X %>% dplyr::select(-participant_ID)
tr_rf3 <- train(Xsub, Y, method='rf',
trControl=trainControl("cv"), importance = TRUE, tuneLength=10)
p3 <- ggplot(varImp(tr_rf3))Put the two figures side by side
grid.arrange(p2,p3, ncol=2)Not surprisingly, if participant_ID is dropped, then
length, weight and age get
assigned higher variable importance.
Our analysis has two objectives:
we wish to know which variables influence the chances of falling;
we wish to predict the threshold angular velocity at a certain speed for new participants.
First, we standardise (mean centered, scaled to have unit standard deviation) the numerical variables. This helps to interpret the magnitude of the estimated coefficients.
preProcValues <- preProcess(d, method = c("center", "scale"))
dTransformed <- predict(preProcValues, d) %>% drop_na()We drop the NA-rows (there are not too many), as further ahead these will be dropped anyway in bas.glm.
Note that not all participants participate in an equal number of experiments:
dTransformed %>% group_by(participant_ID) %>% count() %>%
ggplot(aes(x=participant_ID,y=n)) + geom_col()We use Bayesian model averaging using the BAS-library to
gain insight into which variables are important to be included. We force
categorical variables to be fully included or not (hence we don’t allow
certain factor levels to be included, while other levels of the same
factors not). We insist on a model that contains
angular_momentum, for else we cannot determine the required
threshold value from the inferred model.
We use a uniform prior over all possible models and use the default mixture of Zellner-g-prior on the coefficients.
IT <- 20000 # nr of MCMC iterations
set.seed(15)
bma0 = bas.glm(outcome ~ ., data=dTransformed, method="MCMC",
MCMC.iterations=IT, family=binomial(),
include.always = ~ angular_momentum,
modelprior=uniform(), force.heredity=TRUE)
print(summary(bma0),digits=2)## P(B != 0 | Y) model 1 model 2 model 3 model 4 model 5
## Intercept 1.00 1.000 1.000 1.000 1.000 1.000
## participant_ID16 1.00 1.000 1.000 1.000 1.000 1.000
## participant_ID19 1.00 1.000 1.000 1.000 1.000 1.000
## participant_ID22 1.00 1.000 1.000 1.000 1.000 1.000
## participant_ID6 1.00 1.000 1.000 1.000 1.000 1.000
## participant_ID5 1.00 1.000 1.000 1.000 1.000 1.000
## participant_ID12 1.00 1.000 1.000 1.000 1.000 1.000
## participant_ID2 1.00 1.000 1.000 1.000 1.000 1.000
## participant_ID15 1.00 1.000 1.000 1.000 1.000 1.000
## participant_ID20 1.00 1.000 1.000 1.000 1.000 1.000
## participant_ID11 1.00 1.000 1.000 1.000 1.000 1.000
## participant_ID21 1.00 1.000 1.000 1.000 1.000 1.000
## participant_ID1 1.00 1.000 1.000 1.000 1.000 1.000
## participant_ID10 1.00 1.000 1.000 1.000 1.000 1.000
## participant_ID23 1.00 1.000 1.000 1.000 1.000 1.000
## participant_ID24 1.00 1.000 1.000 1.000 1.000 1.000
## participant_ID13 1.00 1.000 1.000 1.000 1.000 1.000
## participant_ID7 1.00 1.000 1.000 1.000 1.000 1.000
## participant_ID14 1.00 1.000 1.000 1.000 1.000 1.000
## participant_ID4 1.00 1.000 1.000 1.000 1.000 1.000
## participant_ID17 1.00 1.000 1.000 1.000 1.000 1.000
## participant_ID8 1.00 1.000 1.000 1.000 1.000 1.000
## participant_ID18 1.00 1.000 1.000 1.000 1.000 1.000
## participant_ID3 1.00 1.000 1.000 1.000 1.000 1.000
## gendermale 0.13 0.000 0.000 0.000 0.000 0.000
## age 0.17 0.000 0.000 0.000 0.000 0.000
## weight 0.13 0.000 0.000 0.000 0.000 0.000
## length 0.15 0.000 0.000 0.000 0.000 0.000
## skill_performance 0.90 1.000 1.000 1.000 1.000 1.000
## skill_effort 0.30 0.000 0.000 0.000 1.000 0.000
## velocity 1.00 1.000 1.000 1.000 1.000 1.000
## pull_ID 1.00 1.000 1.000 1.000 1.000 1.000
## pull_directionCW 0.67 1.000 0.000 1.000 1.000 1.000
## angular_momentum 1.00 1.000 1.000 1.000 1.000 1.000
## reaction_time 0.20 0.000 0.000 0.000 0.000 0.000
## phi_zero 0.55 1.000 1.000 0.000 1.000 0.000
## phi_dot_zero 0.20 0.000 0.000 0.000 0.000 0.000
## delta_zero 0.29 0.000 0.000 1.000 0.000 0.000
## delta_dot_zero 0.95 1.000 1.000 1.000 1.000 1.000
## psi_zero 0.18 0.000 0.000 0.000 0.000 0.000
## Q_zero 0.26 0.000 0.000 0.000 0.000 0.000
## BF NA 1.000 0.748 0.368 0.324 0.492
## PostProbs NA 0.031 0.027 0.023 0.017 0.016
## R2 NA NA NA NA NA NA
## dim NA 31.000 30.000 31.000 32.000 30.000
## logmarg NA -317.066 -317.356 -318.066 -318.192 -317.776
diagnostics(bma0) # should be on straight lineExperiments containing missing values are dropped in the analysis. As seen in the exploratory data analysis, the fraction of missing data is modest. The following figure visualises the models which get highest posterior probability. Black blocks denote that a variable is not contained in the model
image(bma0, rotate=F)One way to choose a model, is by taking the model for which the marginal posterior inclusion probabilities exceed 1/2, the so called posterior median model.
aa <- tibble(probne0=bma0$probne0, names=bma0$namesx)
partic_p <- aa %>% filter(names=="participant_ID1") %>% select(probne0)
p_inclusion <- aa %>%
filter(!str_detect(names,"participant_ID")) %>%
filter(!str_detect(names,"Interc")) %>%
add_row(names="participant_ID", probne0=as.numeric(partic_p)) %>%
mutate(names2= reorder(names, probne0, median)) %>%
ggplot(aes(x=probne0, y=names2)) +
geom_point() + geom_vline(xintercept=0.5, colour="red") +
labs(x="posterior inclusion probability", y="")
p_inclusionThis suggests a model containing: angular_momentum,
participant_ID, delta_ dot_ zero,
pull_ID, velocity,
skill_performance, pull_direction and
phi_zero.
Now suppose we drop participant_ID from the list of
predictors.
set.seed(13)
bma1 = bas.glm(outcome ~ ., data=dTransformed[,-1], method="MCMC",
include.always = ~ angular_momentum,
MCMC.iterations=IT, family=binomial(),
modelprior=uniform(), force.heredity=TRUE)
print(summary(bma1),digits=2)## P(B != 0 | Y) model 1 model 2 model 3 model 4 model 5
## Intercept 1.00 1.000 1.000 1.00 1.000 1.000
## gendermale 0.99 1.000 1.000 1.00 1.000 1.000
## age 1.00 1.000 1.000 1.00 1.000 1.000
## weight 0.40 0.000 0.000 0.00 0.000 0.000
## length 0.56 1.000 1.000 1.00 1.000 0.000
## skill_performance 0.40 0.000 1.000 0.00 1.000 0.000
## skill_effort 0.20 0.000 0.000 0.00 0.000 0.000
## velocity 1.00 1.000 1.000 1.00 1.000 1.000
## pull_ID 0.98 1.000 1.000 1.00 1.000 1.000
## pull_directionCW 0.54 0.000 0.000 1.00 1.000 0.000
## angular_momentum 1.00 1.000 1.000 1.00 1.000 1.000
## reaction_time 0.13 0.000 0.000 0.00 0.000 0.000
## phi_zero 0.39 0.000 0.000 1.00 0.000 0.000
## phi_dot_zero 0.57 0.000 0.000 1.00 1.000 0.000
## delta_zero 0.91 1.000 1.000 1.00 1.000 1.000
## delta_dot_zero 0.91 1.000 1.000 1.00 1.000 1.000
## psi_zero 0.12 0.000 0.000 0.00 0.000 0.000
## Q_zero 0.93 1.000 1.000 1.00 1.000 1.000
## BF NA 0.851 0.875 1.00 0.887 0.830
## PostProbs NA 0.027 0.021 0.02 0.019 0.018
## R2 NA NA NA NA NA NA
## dim NA 10.000 11.000 13.00 13.000 9.000
## logmarg NA -419.657 -419.629 -419.49 -419.615 -419.682
diagnostics(bma1) # should be on straight line image(bma1, rotate=F)aa1 <- tibble(probne0=bma1$probne0, names=bma1$namesx)
partic_p <- aa1 %>% filter(names=="participant_ID1") %>% select(probne0)
p_inclusion1 <- aa1 %>%
filter(!str_detect(names,"participant_ID")) %>%
filter(!str_detect(names,"Interc")) %>%
add_row(names="participant_ID", probne0=as.numeric(partic_p)) %>%
mutate(names2= reorder(names, probne0, median)) %>%
ggplot(aes(x=probne0, y=names2)) +
geom_point() + geom_vline(xintercept=0.5, colour="red") +
labs(x="posterior inclusion probability", y="")
p_inclusion1Not surprisingly, in this case some participant characteristics enter the median probability model.
We make a plot of confidence intervals for the coefficients, excluding those for participants.
c0 <- confint(coef(bma0))
c0_df <- tibble(coefficient=names(c0[,1]), mean=as.vector(c0[,3]),
low=as.vector(c0[,1]), up= as.vector(c0[,2])) %>%
mutate(type=str_detect(coefficient,"partic"))
c1 <- confint(coef(bma1))
c1_df <- tibble(coefficient=names(c1[,1]), mean=as.vector(c1[,3]),
low=as.vector(c1[,1]), up= as.vector(c1[,2])) %>%
mutate(type=str_detect(coefficient,"partic"))
c01_df <- rbind(c0_df, c1_df) %>% filter(type==FALSE) %>%
mutate(participant_included=rep(c("yes", "no"),each=18))
c01_df %>%
ggplot(aes(x=coefficient, y=mean,color=participant_included)) +
geom_point() +
geom_errorbar(aes(ymin=low, ymax=up)) +
geom_hline(yintercept = 0)+ coord_flip() +
theme(legend.position = "bottom")The preceding analysis has been helpful in identifying important
predictors. It does however neglect the multilevel structure. It does
not seem easy to specify that in bas.glm. For that reason, we refit the
model, including the multilevel structure, using the brms
package.
We first set the prior distributions on the coefficients to be zero-mean Normally distributed with standard-deviation of 2.
prior1 <- set_prior("normal(0, 2)", class = "b")We fit a multilevel version based on the median probability model
appearing from bma0.
set.seed(16)
fit0 <- brm(outcome ~ 0 + Intercept + phi_zero + pull_direction +
skill_performance + delta_dot_zero + velocity +
pull_ID + angular_momentum + (1 | participant_ID),
data=dTransformed,
family = bernoulli(link="logit"),
prior=prior1,
warmup = 5000,
iter = IT,
chains = 4,
seed = 123)## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
##
## SAMPLING FOR MODEL '1919117ad03e75562c036424e5d8ac5f' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.000167 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.67 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 20000 [ 0%] (Warmup)
## Chain 1: Iteration: 2000 / 20000 [ 10%] (Warmup)
## Chain 1: Iteration: 4000 / 20000 [ 20%] (Warmup)
## Chain 1: Iteration: 5001 / 20000 [ 25%] (Sampling)
## Chain 1: Iteration: 7000 / 20000 [ 35%] (Sampling)
## Chain 1: Iteration: 9000 / 20000 [ 45%] (Sampling)
## Chain 1: Iteration: 11000 / 20000 [ 55%] (Sampling)
## Chain 1: Iteration: 13000 / 20000 [ 65%] (Sampling)
## Chain 1: Iteration: 15000 / 20000 [ 75%] (Sampling)
## Chain 1: Iteration: 17000 / 20000 [ 85%] (Sampling)
## Chain 1: Iteration: 19000 / 20000 [ 95%] (Sampling)
## Chain 1: Iteration: 20000 / 20000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 8.15245 seconds (Warm-up)
## Chain 1: 31.6577 seconds (Sampling)
## Chain 1: 39.8101 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL '1919117ad03e75562c036424e5d8ac5f' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 7.6e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.76 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 20000 [ 0%] (Warmup)
## Chain 2: Iteration: 2000 / 20000 [ 10%] (Warmup)
## Chain 2: Iteration: 4000 / 20000 [ 20%] (Warmup)
## Chain 2: Iteration: 5001 / 20000 [ 25%] (Sampling)
## Chain 2: Iteration: 7000 / 20000 [ 35%] (Sampling)
## Chain 2: Iteration: 9000 / 20000 [ 45%] (Sampling)
## Chain 2: Iteration: 11000 / 20000 [ 55%] (Sampling)
## Chain 2: Iteration: 13000 / 20000 [ 65%] (Sampling)
## Chain 2: Iteration: 15000 / 20000 [ 75%] (Sampling)
## Chain 2: Iteration: 17000 / 20000 [ 85%] (Sampling)
## Chain 2: Iteration: 19000 / 20000 [ 95%] (Sampling)
## Chain 2: Iteration: 20000 / 20000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 7.94158 seconds (Warm-up)
## Chain 2: 32.8181 seconds (Sampling)
## Chain 2: 40.7597 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL '1919117ad03e75562c036424e5d8ac5f' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 6.5e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.65 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 20000 [ 0%] (Warmup)
## Chain 3: Iteration: 2000 / 20000 [ 10%] (Warmup)
## Chain 3: Iteration: 4000 / 20000 [ 20%] (Warmup)
## Chain 3: Iteration: 5001 / 20000 [ 25%] (Sampling)
## Chain 3: Iteration: 7000 / 20000 [ 35%] (Sampling)
## Chain 3: Iteration: 9000 / 20000 [ 45%] (Sampling)
## Chain 3: Iteration: 11000 / 20000 [ 55%] (Sampling)
## Chain 3: Iteration: 13000 / 20000 [ 65%] (Sampling)
## Chain 3: Iteration: 15000 / 20000 [ 75%] (Sampling)
## Chain 3: Iteration: 17000 / 20000 [ 85%] (Sampling)
## Chain 3: Iteration: 19000 / 20000 [ 95%] (Sampling)
## Chain 3: Iteration: 20000 / 20000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 7.77018 seconds (Warm-up)
## Chain 3: 30.4534 seconds (Sampling)
## Chain 3: 38.2236 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL '1919117ad03e75562c036424e5d8ac5f' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 9e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.9 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 20000 [ 0%] (Warmup)
## Chain 4: Iteration: 2000 / 20000 [ 10%] (Warmup)
## Chain 4: Iteration: 4000 / 20000 [ 20%] (Warmup)
## Chain 4: Iteration: 5001 / 20000 [ 25%] (Sampling)
## Chain 4: Iteration: 7000 / 20000 [ 35%] (Sampling)
## Chain 4: Iteration: 9000 / 20000 [ 45%] (Sampling)
## Chain 4: Iteration: 11000 / 20000 [ 55%] (Sampling)
## Chain 4: Iteration: 13000 / 20000 [ 65%] (Sampling)
## Chain 4: Iteration: 15000 / 20000 [ 75%] (Sampling)
## Chain 4: Iteration: 17000 / 20000 [ 85%] (Sampling)
## Chain 4: Iteration: 19000 / 20000 [ 95%] (Sampling)
## Chain 4: Iteration: 20000 / 20000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 8.56809 seconds (Warm-up)
## Chain 4: 30.9096 seconds (Sampling)
## Chain 4: 39.4777 seconds (Total)
## Chain 4:
fit0 %>% summary()## Family: bernoulli
## Links: mu = logit
## Formula: outcome ~ 0 + Intercept + phi_zero + pull_direction + skill_performance + delta_dot_zero + velocity + pull_ID + angular_momentum + (1 | participant_ID)
## Data: dTransformed (Number of observations: 874)
## Draws: 4 chains, each with iter = 20000; warmup = 5000; thin = 1;
## total post-warmup draws = 60000
##
## Group-Level Effects:
## ~participant_ID (Number of levels: 24)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 3.69 0.62 2.67 5.07 1.00 11133 20445
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.48 0.72 -0.93 1.94 1.00 9167 15700
## phi_zero 0.42 0.18 0.07 0.77 1.00 41696 41662
## pull_directionCW -0.49 0.25 -0.98 -0.01 1.00 47048 42765
## skill_performance 0.71 0.22 0.28 1.14 1.00 31549 39407
## delta_dot_zero -0.49 0.15 -0.79 -0.18 1.00 43822 42088
## velocity 1.62 0.18 1.27 1.99 1.00 33574 39211
## pull_ID -0.77 0.15 -1.06 -0.48 1.00 43260 40636
## angular_momentum 5.26 0.38 4.53 6.03 1.00 33461 36826
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
prior_summary(fit0)## prior class coef group resp dpar nlpar
## normal(0, 2) b
## normal(0, 2) b angular_momentum
## normal(0, 2) b delta_dot_zero
## normal(0, 2) b Intercept
## normal(0, 2) b phi_zero
## normal(0, 2) b pull_directionCW
## normal(0, 2) b pull_ID
## normal(0, 2) b skill_performance
## normal(0, 2) b velocity
## student_t(3, 0, 2.5) sd
## student_t(3, 0, 2.5) sd participant_ID
## student_t(3, 0, 2.5) sd Intercept participant_ID
## bound source
## user
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## default
## (vectorized)
## (vectorized)
These coefficients are not exactly the same as under bma0 since
we don’t use Bayesian model averaging;
the prior distribution is different;
we employ a multilevel model.
dTransformed_reformat <- dTransformed %>%
mutate(outcomenum = as.numeric(outcome)-1)
glmer_fit0 <- glmer(outcomenum ~ phi_zero + pull_direction +
skill_performance + delta_dot_zero + velocity +
pull_ID + angular_momentum + (1 | participant_ID),
data=dTransformed_reformat, family = binomial)
glmer_fit0 %>% summary()## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: outcomenum ~ phi_zero + pull_direction + skill_performance +
## delta_dot_zero + velocity + pull_ID + angular_momentum +
## (1 | participant_ID)
## Data: dTransformed_reformat
##
## AIC BIC logLik deviance df.resid
## 624.1 667.0 -303.0 606.1 865
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.6102 -0.2992 -0.0160 0.2774 7.3427
##
## Random effects:
## Groups Name Variance Std.Dev.
## participant_ID (Intercept) 13.28 3.644
## Number of obs: 874, groups: participant_ID, 24
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.5533 0.7689 0.720 0.47180
## phi_zero 0.4236 0.1795 2.360 0.01828 *
## pull_directionCW -0.5070 0.2500 -2.028 0.04257 *
## skill_performance 0.7343 0.2252 3.261 0.00111 **
## delta_dot_zero -0.4956 0.1561 -3.174 0.00150 **
## velocity 1.6616 0.1879 8.844 < 2e-16 ***
## pull_ID -0.7837 0.1527 -5.132 2.87e-07 ***
## angular_momentum 5.3777 0.4015 13.394 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) phi_zr pll_CW skll_p dlt_d_ velcty pll_ID
## phi_zero -0.033
## pll_drctnCW -0.177 0.078
## skll_prfrmn -0.023 0.049 -0.030
## delta_dt_zr 0.059 -0.258 -0.388 -0.014
## velocity 0.032 0.084 -0.098 0.477 -0.060
## pull_ID 0.010 -0.033 0.023 0.267 0.125 -0.024
## anglr_mmntm 0.049 0.120 -0.099 0.240 -0.205 0.608 -0.397
To facilitate comparison, we show both the Bayesian and frequentist fitted models:
sum0 <- fit0 %>% summary() %>% unclass()
sum0$fixed## Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept 0.4829009 0.7246277 -0.93089505 1.93835361 1.000372
## phi_zero 0.4169201 0.1770367 0.07206348 0.76711668 1.000102
## pull_directionCW -0.4937041 0.2457098 -0.97539265 -0.01218042 1.000055
## skill_performance 0.7059425 0.2203738 0.27837368 1.14201295 1.000031
## delta_dot_zero -0.4865379 0.1541700 -0.78841881 -0.18386285 1.000090
## velocity 1.6197151 0.1828174 1.27426927 1.99044089 1.000170
## pull_ID -0.7658691 0.1495807 -1.06476398 -0.47598085 1.000059
## angular_momentum 5.2557019 0.3801358 4.53472387 6.02727389 1.000292
## Bulk_ESS Tail_ESS
## Intercept 9166.944 15699.62
## phi_zero 41696.270 41662.04
## pull_directionCW 47047.902 42765.20
## skill_performance 31549.339 39407.11
## delta_dot_zero 43821.940 42087.76
## velocity 33574.146 39211.02
## pull_ID 43259.604 40636.39
## angular_momentum 33460.657 36825.72
sum0$random## $participant_ID
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 3.694535 0.6165921 2.673694 5.07287 1.0003 11133.41 20445.24
to be compared to
glmer0 <- glmer_fit0 %>% summary() %>% unclass()
coefficients(glmer0)## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.5532895 0.7689346 0.7195534 4.718000e-01
## phi_zero 0.4236194 0.1795149 2.3598010 1.828474e-02
## pull_directionCW -0.5070225 0.2500183 -2.0279420 4.256617e-02
## skill_performance 0.7342650 0.2251922 3.2606145 1.111711e-03
## delta_dot_zero -0.4955544 0.1561312 -3.1739608 1.503740e-03
## velocity 1.6616009 0.1878872 8.8436099 9.267472e-19
## pull_ID -0.7837200 0.1527248 -5.1315842 2.873137e-07
## angular_momentum 5.3777007 0.4015094 13.3937119 6.580788e-41
glmer0$varcor## Groups Name Std.Dev.
## participant_ID (Intercept) 3.644
The posterior means from brm and estimates from
glmer are very similar. Additionally, we obtain
significance of all variables in the model.
Posterior predictive checks are used to establish visually if the fitted model captures certain features in the data appropriately. Here we define such a check.
The function checks2 compares the fraction of falls for
the total number of perturbations of each participant to that based on
1000 predictions of the fitted model. It does so separately for each of
the velocities, which are in the set {6,12,18}
checks2 <- function(d, fit) # fraction falling by participant and factor(velocity)
{
pred <- posterior_predict(fit, ndraws=1000)
fracfalling1 <- d %>%
mutate(velocity=as.factor(round(velocity,1))) %>%
mutate(predicted_average=colMeans(pred)) %>%
group_by(velocity, participant_ID) %>%
summarise(n=n(), observed=mean(outcome=="1"),
predicted=mean(predicted_average)) %>%
pivot_longer(cols=c(predicted, observed), names_to="type", values_to="fraction")
p <- fracfalling1 %>%
ggplot() +
geom_point(aes(x=velocity, y = fraction,colour=type)) +
labs(x="velocity (standardized)", y="fraction of falls of the total number of perturbations")
p + facet_wrap(~participant_ID)
}
checks2(dTransformed, fit0)We check the effect of the specified prior distribution on the estimates for the coefficients.
prior_vague <- set_prior("normal(0, 5)", class = "b") # way larger sd
prior_centred <- set_prior("normal(0, 1)", class = "b")
# somewhat small sd, possibly pulls estimates towards zero.Also fit the model with these priors:
fit0_vague <- brm(outcome ~ 0 + Intercept + phi_zero + pull_direction +
skill_performance + delta_dot_zero + velocity +
pull_ID + angular_momentum + (1 | participant_ID),
data=dTransformed,
family = bernoulli(link="logit"),
prior=prior_vague,
warmup = 5000,
iter = IT,
chains = 4,
seed = 123)## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
##
## SAMPLING FOR MODEL 'ce805167e8faadb918648a9d001dabfe' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.000224 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.24 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 20000 [ 0%] (Warmup)
## Chain 1: Iteration: 2000 / 20000 [ 10%] (Warmup)
## Chain 1: Iteration: 4000 / 20000 [ 20%] (Warmup)
## Chain 1: Iteration: 5001 / 20000 [ 25%] (Sampling)
## Chain 1: Iteration: 7000 / 20000 [ 35%] (Sampling)
## Chain 1: Iteration: 9000 / 20000 [ 45%] (Sampling)
## Chain 1: Iteration: 11000 / 20000 [ 55%] (Sampling)
## Chain 1: Iteration: 13000 / 20000 [ 65%] (Sampling)
## Chain 1: Iteration: 15000 / 20000 [ 75%] (Sampling)
## Chain 1: Iteration: 17000 / 20000 [ 85%] (Sampling)
## Chain 1: Iteration: 19000 / 20000 [ 95%] (Sampling)
## Chain 1: Iteration: 20000 / 20000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 8.72305 seconds (Warm-up)
## Chain 1: 32.1648 seconds (Sampling)
## Chain 1: 40.8878 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'ce805167e8faadb918648a9d001dabfe' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 8.3e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.83 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 20000 [ 0%] (Warmup)
## Chain 2: Iteration: 2000 / 20000 [ 10%] (Warmup)
## Chain 2: Iteration: 4000 / 20000 [ 20%] (Warmup)
## Chain 2: Iteration: 5001 / 20000 [ 25%] (Sampling)
## Chain 2: Iteration: 7000 / 20000 [ 35%] (Sampling)
## Chain 2: Iteration: 9000 / 20000 [ 45%] (Sampling)
## Chain 2: Iteration: 11000 / 20000 [ 55%] (Sampling)
## Chain 2: Iteration: 13000 / 20000 [ 65%] (Sampling)
## Chain 2: Iteration: 15000 / 20000 [ 75%] (Sampling)
## Chain 2: Iteration: 17000 / 20000 [ 85%] (Sampling)
## Chain 2: Iteration: 19000 / 20000 [ 95%] (Sampling)
## Chain 2: Iteration: 20000 / 20000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 8.37378 seconds (Warm-up)
## Chain 2: 32.1127 seconds (Sampling)
## Chain 2: 40.4865 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'ce805167e8faadb918648a9d001dabfe' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 5.8e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.58 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 20000 [ 0%] (Warmup)
## Chain 3: Iteration: 2000 / 20000 [ 10%] (Warmup)
## Chain 3: Iteration: 4000 / 20000 [ 20%] (Warmup)
## Chain 3: Iteration: 5001 / 20000 [ 25%] (Sampling)
## Chain 3: Iteration: 7000 / 20000 [ 35%] (Sampling)
## Chain 3: Iteration: 9000 / 20000 [ 45%] (Sampling)
## Chain 3: Iteration: 11000 / 20000 [ 55%] (Sampling)
## Chain 3: Iteration: 13000 / 20000 [ 65%] (Sampling)
## Chain 3: Iteration: 15000 / 20000 [ 75%] (Sampling)
## Chain 3: Iteration: 17000 / 20000 [ 85%] (Sampling)
## Chain 3: Iteration: 19000 / 20000 [ 95%] (Sampling)
## Chain 3: Iteration: 20000 / 20000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 9.01871 seconds (Warm-up)
## Chain 3: 33.0101 seconds (Sampling)
## Chain 3: 42.0288 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'ce805167e8faadb918648a9d001dabfe' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 6e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.6 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 20000 [ 0%] (Warmup)
## Chain 4: Iteration: 2000 / 20000 [ 10%] (Warmup)
## Chain 4: Iteration: 4000 / 20000 [ 20%] (Warmup)
## Chain 4: Iteration: 5001 / 20000 [ 25%] (Sampling)
## Chain 4: Iteration: 7000 / 20000 [ 35%] (Sampling)
## Chain 4: Iteration: 9000 / 20000 [ 45%] (Sampling)
## Chain 4: Iteration: 11000 / 20000 [ 55%] (Sampling)
## Chain 4: Iteration: 13000 / 20000 [ 65%] (Sampling)
## Chain 4: Iteration: 15000 / 20000 [ 75%] (Sampling)
## Chain 4: Iteration: 17000 / 20000 [ 85%] (Sampling)
## Chain 4: Iteration: 19000 / 20000 [ 95%] (Sampling)
## Chain 4: Iteration: 20000 / 20000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 8.67467 seconds (Warm-up)
## Chain 4: 31.6836 seconds (Sampling)
## Chain 4: 40.3583 seconds (Total)
## Chain 4:
fit0_centred <- brm(outcome ~ 0 + Intercept + phi_zero + pull_direction +
skill_performance + delta_dot_zero + velocity +
pull_ID + angular_momentum + (1 | participant_ID),
data=dTransformed,
family = bernoulli(link="logit"),
prior=prior_centred,
warmup = 5000,
iter = IT,
chains = 4,
seed = 123)## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
##
## SAMPLING FOR MODEL 'a07710f062e3c373f1c44960e5a29b2c' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.000155 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.55 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 20000 [ 0%] (Warmup)
## Chain 1: Iteration: 2000 / 20000 [ 10%] (Warmup)
## Chain 1: Iteration: 4000 / 20000 [ 20%] (Warmup)
## Chain 1: Iteration: 5001 / 20000 [ 25%] (Sampling)
## Chain 1: Iteration: 7000 / 20000 [ 35%] (Sampling)
## Chain 1: Iteration: 9000 / 20000 [ 45%] (Sampling)
## Chain 1: Iteration: 11000 / 20000 [ 55%] (Sampling)
## Chain 1: Iteration: 13000 / 20000 [ 65%] (Sampling)
## Chain 1: Iteration: 15000 / 20000 [ 75%] (Sampling)
## Chain 1: Iteration: 17000 / 20000 [ 85%] (Sampling)
## Chain 1: Iteration: 19000 / 20000 [ 95%] (Sampling)
## Chain 1: Iteration: 20000 / 20000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 7.19273 seconds (Warm-up)
## Chain 1: 26.2807 seconds (Sampling)
## Chain 1: 33.4734 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'a07710f062e3c373f1c44960e5a29b2c' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 6.9e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.69 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 20000 [ 0%] (Warmup)
## Chain 2: Iteration: 2000 / 20000 [ 10%] (Warmup)
## Chain 2: Iteration: 4000 / 20000 [ 20%] (Warmup)
## Chain 2: Iteration: 5001 / 20000 [ 25%] (Sampling)
## Chain 2: Iteration: 7000 / 20000 [ 35%] (Sampling)
## Chain 2: Iteration: 9000 / 20000 [ 45%] (Sampling)
## Chain 2: Iteration: 11000 / 20000 [ 55%] (Sampling)
## Chain 2: Iteration: 13000 / 20000 [ 65%] (Sampling)
## Chain 2: Iteration: 15000 / 20000 [ 75%] (Sampling)
## Chain 2: Iteration: 17000 / 20000 [ 85%] (Sampling)
## Chain 2: Iteration: 19000 / 20000 [ 95%] (Sampling)
## Chain 2: Iteration: 20000 / 20000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 6.8328 seconds (Warm-up)
## Chain 2: 24.6561 seconds (Sampling)
## Chain 2: 31.4889 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'a07710f062e3c373f1c44960e5a29b2c' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 5.9e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.59 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 20000 [ 0%] (Warmup)
## Chain 3: Iteration: 2000 / 20000 [ 10%] (Warmup)
## Chain 3: Iteration: 4000 / 20000 [ 20%] (Warmup)
## Chain 3: Iteration: 5001 / 20000 [ 25%] (Sampling)
## Chain 3: Iteration: 7000 / 20000 [ 35%] (Sampling)
## Chain 3: Iteration: 9000 / 20000 [ 45%] (Sampling)
## Chain 3: Iteration: 11000 / 20000 [ 55%] (Sampling)
## Chain 3: Iteration: 13000 / 20000 [ 65%] (Sampling)
## Chain 3: Iteration: 15000 / 20000 [ 75%] (Sampling)
## Chain 3: Iteration: 17000 / 20000 [ 85%] (Sampling)
## Chain 3: Iteration: 19000 / 20000 [ 95%] (Sampling)
## Chain 3: Iteration: 20000 / 20000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 7.05009 seconds (Warm-up)
## Chain 3: 23.7962 seconds (Sampling)
## Chain 3: 30.8463 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'a07710f062e3c373f1c44960e5a29b2c' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 6.1e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.61 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 20000 [ 0%] (Warmup)
## Chain 4: Iteration: 2000 / 20000 [ 10%] (Warmup)
## Chain 4: Iteration: 4000 / 20000 [ 20%] (Warmup)
## Chain 4: Iteration: 5001 / 20000 [ 25%] (Sampling)
## Chain 4: Iteration: 7000 / 20000 [ 35%] (Sampling)
## Chain 4: Iteration: 9000 / 20000 [ 45%] (Sampling)
## Chain 4: Iteration: 11000 / 20000 [ 55%] (Sampling)
## Chain 4: Iteration: 13000 / 20000 [ 65%] (Sampling)
## Chain 4: Iteration: 15000 / 20000 [ 75%] (Sampling)
## Chain 4: Iteration: 17000 / 20000 [ 85%] (Sampling)
## Chain 4: Iteration: 19000 / 20000 [ 95%] (Sampling)
## Chain 4: Iteration: 20000 / 20000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 7.07604 seconds (Warm-up)
## Chain 4: 22.6028 seconds (Sampling)
## Chain 4: 29.6788 seconds (Total)
## Chain 4:
fit0 %>% summary()## Family: bernoulli
## Links: mu = logit
## Formula: outcome ~ 0 + Intercept + phi_zero + pull_direction + skill_performance + delta_dot_zero + velocity + pull_ID + angular_momentum + (1 | participant_ID)
## Data: dTransformed (Number of observations: 874)
## Draws: 4 chains, each with iter = 20000; warmup = 5000; thin = 1;
## total post-warmup draws = 60000
##
## Group-Level Effects:
## ~participant_ID (Number of levels: 24)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 3.69 0.62 2.67 5.07 1.00 11133 20445
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.48 0.72 -0.93 1.94 1.00 9167 15700
## phi_zero 0.42 0.18 0.07 0.77 1.00 41696 41662
## pull_directionCW -0.49 0.25 -0.98 -0.01 1.00 47048 42765
## skill_performance 0.71 0.22 0.28 1.14 1.00 31549 39407
## delta_dot_zero -0.49 0.15 -0.79 -0.18 1.00 43822 42088
## velocity 1.62 0.18 1.27 1.99 1.00 33574 39211
## pull_ID -0.77 0.15 -1.06 -0.48 1.00 43260 40636
## angular_momentum 5.26 0.38 4.53 6.03 1.00 33461 36826
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
fit0_vague %>% summary()## Family: bernoulli
## Links: mu = logit
## Formula: outcome ~ 0 + Intercept + phi_zero + pull_direction + skill_performance + delta_dot_zero + velocity + pull_ID + angular_momentum + (1 | participant_ID)
## Data: dTransformed (Number of observations: 874)
## Draws: 4 chains, each with iter = 20000; warmup = 5000; thin = 1;
## total post-warmup draws = 60000
##
## Group-Level Effects:
## ~participant_ID (Number of levels: 24)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 3.84 0.65 2.78 5.30 1.00 11541 20586
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.57 0.80 -0.99 2.14 1.00 9069 15652
## phi_zero 0.43 0.18 0.08 0.79 1.00 41898 40680
## pull_directionCW -0.51 0.25 -1.01 -0.03 1.00 46950 40973
## skill_performance 0.74 0.23 0.31 1.19 1.00 31771 37675
## delta_dot_zero -0.50 0.16 -0.81 -0.20 1.00 44348 42126
## velocity 1.68 0.19 1.33 2.06 1.00 32655 36795
## pull_ID -0.79 0.15 -1.10 -0.50 1.00 41567 39251
## angular_momentum 5.46 0.40 4.71 6.28 1.00 31567 34596
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
fit0_centred %>% summary()## Family: bernoulli
## Links: mu = logit
## Formula: outcome ~ 0 + Intercept + phi_zero + pull_direction + skill_performance + delta_dot_zero + velocity + pull_ID + angular_momentum + (1 | participant_ID)
## Data: dTransformed (Number of observations: 874)
## Draws: 4 chains, each with iter = 20000; warmup = 5000; thin = 1;
## total post-warmup draws = 60000
##
## Group-Level Effects:
## ~participant_ID (Number of levels: 24)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 3.32 0.55 2.41 4.57 1.00 11275 22038
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.33 0.58 -0.82 1.45 1.00 8446 15627
## phi_zero 0.38 0.17 0.05 0.71 1.00 36554 40778
## pull_directionCW -0.44 0.23 -0.89 0.01 1.00 44925 42130
## skill_performance 0.61 0.20 0.22 1.02 1.00 29294 37257
## delta_dot_zero -0.44 0.15 -0.73 -0.16 1.00 43493 44123
## velocity 1.45 0.17 1.13 1.78 1.00 32566 37339
## pull_ID -0.69 0.14 -0.97 -0.42 1.00 41114 39827
## angular_momentum 4.73 0.32 4.13 5.38 1.00 31855 39929
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
From here, we see that the differences with the “vague” prior are small; the “centred” prior shrinks coefficients more towards zero and appears to have some effect on the estimates.
Some example plots to illustrate the effect of a particular variable
on outcome.
conditional_effects(fit0, effects = "pull_direction")conditional_effects(fit0, effects = c("velocity"))conditional_effects(fit0, effects = c("delta_dot_zero"))conditional_effects(fit0, effects = c("angular_momentum:velocity"))We visualise the posterior by showing density- and traceplots of posterior samples.
fit0 %>% plot(newpage=FALSE, ask=FALSE)We extract IT posterior samples from the fitted
model.
ps <- brms::as_draws_df(fit0) %>% tibble() %>% sample_n(IT)We plot all parameters, except for participant effects
ln <- c("b_Intercept", "b_skill_performance", "b_velocity", "b_phi_zero",
"b_pull_ID", "b_pull_directionCW", "b_angular_momentum",
"b_delta_dot_zero","sd_participant_ID__Intercept")
ps_coef <- ps %>% dplyr::select(-contains("r_particip")) %>%
dplyr::select(contains("b_") | contains("sd_")) %>%
tidyr::pivot_longer(cols=ln,names_to="coefficient", values_to="value") %>%
mutate(coefficient=str_remove(coefficient,"b_")) %>%
mutate(coefficient = recode(coefficient, `sd_participant_ID__Intercept` = "sd_participant"))
ps_coef %>% ggplot(aes(x=value, y = reorder(coefficient, value,median))) +
stat_halfeye( point_size=1, interval_size=0.5)+
labs(x="coefficient", y = "")Similarly for participant effects
ps_sub <- ps %>% dplyr::select(contains("r_particip"))
ps_sublong <- ps_sub %>% tidyr::pivot_longer(cols=contains("r_particip"),
names_to="participant", values_to="value") %>%
mutate(participant_ID=as.factor(readr::parse_number(participant)))
ps_complete <- left_join(ps_sublong, d, by ="participant_ID")
ps_complete %>% ggplot(aes(x=value,
y = reorder(participant_ID, value,median),fill=age )) +
stat_halfeye(point_colour="orange", point_size=1, interval_size=0.5)+
# labs(x="coefficient", y = "participant_ID")
labs(x=TeX(r'($\beta_i$)'), y="participant_ID")rm(ps_coef, ps_sub, ps_sublong, ps_complete) # free memoryConsider the results for following persons:
participant 1 with 12 km/h and 6 km/h
participant 9 with 12 km/h and 6 km/h
participant 11 with 12 km/h and with 6 km/h and with 18 km/h
participant 12 with 12 km/h and with 6 km/h
phi_zero = 0 (bicycle in straight position),
pull_directionCW = 1, delta_dot_zero = 0,
pull_id = 1
pIDs <- c(1,1,9,9,11,11,11,12,12) # arbitrarily choose a few participantsd1 <- d %>% filter(participant_ID==as.character(1)) %>% slice(1)
d9 <- d %>% filter(participant_ID==as.character(9)) %>% slice(1)
d11 <- d %>% filter(participant_ID==as.character(11)) %>% slice(1)
d12 <- d %>% filter(participant_ID==as.character(12)) %>% slice(1)
d_existing <- bind_rows(d1,d1,d9,d9,d11,d11,d11,d12,d12) %>%
mutate(participant_ID2=c("P1_v6","P1_v12","P9_v6", "P9_v12", "P11_v6","P11_v12","P11_v18", "P12_v6", "P12_v12" )) %>%
mutate(velocity=c(6,12,6,12,6,12,18,6,12)) %>%
mutate(phi_zero=rep(0,9), pull_direction=rep("CW",9), delta_dot_zero=rep(0,9), pull_ID=rep(1,9))
d_existingTransformed <- predict(preProcValues, d_existing)
nps <- nrow(ps) # nr of posterior samples
zz = rep(0,nps)
z <- data.frame("P1_v6" = zz,"P1_v12" = zz,"P9_v6" = zz, "P9_v12" = zz, "P11_v6" = zz,
"P11_v12" = zz,"P11_v18" = zz, "P12_v6" = zz, "P12_v12" = zz)
head(z)## P1_v6 P1_v12 P9_v6 P9_v12 P11_v6 P11_v12 P11_v18 P12_v6 P12_v12
## 1 0 0 0 0 0 0 0 0 0
## 2 0 0 0 0 0 0 0 0 0
## 3 0 0 0 0 0 0 0 0 0
## 4 0 0 0 0 0 0 0 0 0
## 5 0 0 0 0 0 0 0 0 0
## 6 0 0 0 0 0 0 0 0 0
angular_threshold <- function(participant_ch, postdraw, participant_postdraw)
{
xx <- participant_postdraw +
postdraw$b_Intercept +
postdraw$b_skill_performance * participant_ch$skill_performance +
postdraw$b_velocity * participant_ch$velocity+
postdraw$b_pull_ID * participant_ch$pull_ID +
postdraw$b_pull_directionCW * (participant_ch$pull_direction=="CW") +
# postdraw$b_phi_zero * participant_ch$phi_zero +
postdraw$b_delta_dot_zero * participant_ch$delta_dot_zero
out <- -xx/(postdraw$b_angular_momentum)
out
}
for (j in 1:9)
{
ID <- pIDs[j]
# select for the participant posterior draws from its intercept:
lab <- paste0("r_participant_ID[",as.character(ID),",Intercept]")
raneff <- ps[,lab] %>% deframe()
# choose participant experimental setting:
participant_ch <- d_existingTransformed[j,]
for (i in 1:nps) # for each posterior sample
{
postdraw <- ps %>% slice(i)
z[i, j] <- angular_threshold(participant_ch, postdraw, raneff[i])
}
}
m_ang <- mean(na.omit(d)$angular_momentum)
sd_ang <- sd(na.omit(d)$angular_momentum)
angular_threshold_df <- apply(z,2, function(x) { m_ang + sd_ang * x } ) # rescale each column
angular_threshold_tidy <- as_tibble(angular_threshold_df) %>%
pivot_longer(cols=1:9, names_to="participant", values_to="eta") angular_threshold_tidy %>%
#ggplot(aes(x=eta, y=reorder(participant, eta, median))) +
ggplot(aes(x=eta, y=participant)) +
stat_halfeye(fill="lightblue") +
labs(x=TeX(r'($\Delta L_{50\%}$)'), y ="") +
ggtitle("Threshold value for angular momentum")We manually construct “new” participants, for which we will compute the perturbation threshold (for the variables in model 5 we provide values, while all remaining columns are just filled with NA-values).
skill performance 40, velocity 12 km/h
skill performance 90, velocity 12 km/h
skill performance 40, velocity 6 km/h
skill performance 90, velocity 6 km/h
skill performance 40, velocity 18 km/h
skill performance 40, velocity 25 km/h (velocity not included in experiment)
dpred <- tibble(participant_ID=c("new_v12_sp40", "new_v12_sp90", "new_v6_sp40", "new_v6_sp90","new_v18_sp40","new_v25_sp40"),
velocity=c(12,12,6,6,18,25), age = rep(NA,6),
skill_performance=c(40,90,40,90,40,40), pull_ID=rep(1,6),
pull_direction=rep("CW",6), weight=rep(NA,6), length=rep(NA,6),
skill_effort = rep(NA,6), initial_search=rep(NA,6), phi_zero=rep(0,6),
phi_dot_zero=rep(NA,6), delta_zero=rep(NA,6), delta_dot_zero=rep(0,6),
psi_zero=rep(NA,6), Q_zero=rep(NA,6), comments=rep(NA,6),
exclude=rep(NA,6), reaction_time=rep(NA,6), pull_force=rep(NA,6),
angular_momentum=rep(NA,6))Now we will predict for these new participants:
dpred_tf <- predict(preProcValues, dpred) # tf=transformed as we did for the model
dpred_tf$delta_dot_zero <- rep(0,6)
nparticipants = nrow(dpred) # nr of participants for which we wish to predict
zz = rep(0,nps)
z <- data.frame(new_participant1 = zz, new_participant2 = zz,
new_participant3 = zz,new_participant4 = zz,
new_participant5 = zz, new_participant6 = zz)
for (j in 1:nparticipants) # for each participant we wish to find threshold value
{
for (i in 1:nps) # for each participant
{
particip = rnorm(1,sd=ps$sd_participant_ID__Intercept[i])
z[i, j] <- angular_threshold(dpred_tf[j,], ps[i,], particip)
}
}
# rescale each column:
angular_threshold_df_new <- apply(z,2, function(x) { m_ang + sd_ang * x})
angular_threshold_tidy_new <- as_tibble(angular_threshold_df_new) %>%
pivot_longer(cols=contains("new"), names_to="participant", values_to="eta") %>%
mutate(participant=fct_recode(participant,
s40v12="new_participant1",
s90v12="new_participant2",
s40v6="new_participant3",
s90v6="new_participant4",
s40v18="new_participant5",
s40v25="new_participant6"))angular_threshold_tidy_new %>%
# ggplot(aes(x=eta, y=reorder(participant, eta, median))) +
ggplot(aes(x=eta, y=participant)) +
stat_halfeye(fill="lightblue") +
labs(x=TeX(r'($\Delta L_{50\%}$)'), y ="") +
ggtitle("Threshold value for angular momentum")We can also make a combined plot for participants which were part of the study and “new” participants:
nr <- nrow(angular_threshold_tidy)
nr_new <- nrow(angular_threshold_tidy_new)
factor_new <- c(rep("no",nr), rep("yes",nr_new))
bind_rows(angular_threshold_tidy, angular_threshold_tidy_new) %>%
mutate(new=factor_new) %>%
ggplot(aes(x=eta, y=reorder(participant, eta, median), fill=new)) +
stat_halfeye(fill="lightblue") + labs(x="eta", y = "participant") This shows that we can fairly accurately find the threshold value for participants within the study, but have large uncertainty about this value for “future” unseen participants. This is somewhat natural, as we don’t know the inherent cycling performance of unseen participants.